arXiv:1501.04320vl [math.AP] 18 Jan 2015 


Some Free Boundary Problems involving 
Nonlocal Diffusion and Aggregation 

Jose Antonio Carrillo* Juan Luis Vazquez^ 
January 20, 2015 


Keywords: nonlocal diffusion, interaction energy, aggregation, obstacle problems 

Abstract 

We report on recent progress in the study of evolution processes involving degenerate 
parabolic equations what may exhibit free boundaries. The equations we have selected 
follow to recent trends in diffusion theory: considering anomalous diffusion with long- 
range effects, which leads to fractional operators or other operators involving kernels with 
large tails; and the combination of diffusion and aggregation effects, leading to delicate 
long-term equilibria whose description is still incipient. 


1 Introduction 

The systematic mathematical study of free boundary problems is quite recent compared to 
the present importance of the field, and goes back to the last half of the past century. In 
this period it has developed in many directions, which combine modeling of quite diverse 
natural phenomena, particularly in diffusion, wave mechanics and elasticity. A main topic 
in this broad landscape has always been the study of evolution processes where the state 
equations governing the different phases are nonlinear parabolic equations, possibly degenerate 
or singular at certain points or for certain values of the unknowns. This is the case of the 
well-known models like Stefan Problem for the common evolution of two immiscible fluids, 
the Porous Medium Equation, and other models to be mentioned below. 

Free boundary models combine the difficulty of the analysis of the nonlinear PDEs with 
the difficulty in locating the separating interface between the phases, or free boundary in the 
absence of a theory of classical solutions that is generally unavailable. 

So in principle we can solve the problem in some generalized sense but neither the state 
variables are guaranteed to be continuous or the free boundaries to be reasonable hyper¬ 
surfaces. Even when the initial and boundary data are smooth, nonlinear processes involving 
free boundaries may develop singularities in finite time both in the solution and the free 
boundary. This is therefore a difficult problem for both analysis and geometry. In the last 50 
years a considerable progress has been achieved and this is reported in different monographs, 
and in its recent trends in other articles of this volume. Several new models have been actively 
pursued in recent years with various levels of mathematical rigor versus practical applicability. 
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In this article we will discuss free boundary problems arising in two new scenarios, nonlocal 
diffusion and aggregation/swarming processes. 

2 Nonlocal diffusion 

The classical theory of diffusion is expressed mathematically by means of the heat equation, 
and more generally by parabolic equations, normally of the linear type; such approach has 
had an enormous success and is now a foundation stone in science and technology. The last 
half of the past century has witnessed intense activity and progress in the theories of nonlinear 
diffusion, examples being the Stefan Problem, the Porous Medium Equation, the p-Laplacian 
equation, the Total Variation Flow, evolution problems of Hele-Shaw type, the Keller-Segel 
chemotaxis system, and many others. Mean curvature flows also belong to this broad category. 
Reaction diffusion has also attracted considerable attention. 

In the last decade there has been a surge of activity focused on the use of so-called 
fractional diffusion operators to replace the standard Laplace operator (and other kinds of 
elliptic operators with variable coefficients), with the aim of further extending the theory 
by taking into account the presence of the long range interactions that occur in a number 
applications. The new operators do not act by pointwise differentiation but by a global 
integration with respect to a singular kernel; in that way the nonlocal character of the process 
is expressed. 

More generally, research on non-local operators has a tradition related to stochastic pro¬ 
cesses, and has witnessed a rapid expansion in the last decade when it has attracted the 
attention of PDE experts who brought new problems and techniques into the field. The area 
poses new challenges to both pure and applied mathematicians, and is now at the interface 
of at least three wide fields: functional analysis, stochastic processes and partial differential 
equations, while providing a new paradigm in scientific modeling. This interaction between 
mathematics and applications is giving rise to new concepts and methods, and is expected to 
produce new challenging mathematical problems for many years to come. 

Fractional operators 

Though there is a wide class of interesting nonlocal operators under scrutiny, both in theory 
and applications, a substantial part of the current work deals with diffusion modeled by the 
so-called fractional Laplacians. We recall that the fractional Laplacian operator is a kind of 
isotropic differentiation operator of order 2s, for some s € (0,1), that can be conveniently 
defined through its Fourier Transform symbol, which is |£| 2s . Thus, if g is a function in the 
Schwartz class in M^, N > 1, we write (—A ) s g = h if 

h(0 = \tf s 9(0, 

so that for s = 1 we recover the standard Laplacian. This definition allows for a wider range 
of parameters s. The interval of interest for fractional diffusion is 0 < s < 1, and for s < 1 
we can also use the integral representation 

(-A)«g(x) = C N , 2 , p.v. f shl dz< (2.1) 

Jr n \ x ~ z \ 

where P.V. stands for principal value and CV i(T is a normalization constant, with precise value 
that is found in the literature. In the limits s —> 0 and s —> 1 it is possible to recover 
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respectively the identity or the standard minus Laplacian, —A, cf. [21]. Remarkably, the 
latter one cannot be represented by a nonlocal formula of the type (2.1). It is also useful to 
recall that the operators (—A) _s , 0 < s < 1, inverse of the former ones, are given by standard 
convolution expressions: 

(-A)- 9 (x) = C N ,_ 2 » f iz, 

Jr n \ x — z \ 

in terms of the usual Riesz potentials. Basic references for these operators are the books by 
Landkof [68] and Stein [79]. A word of caution: in the literature we often find the notation 
a = 2s, and then the desired interval is 0 < cr < 2. According to that practice, we will 
sometimes use a instead of s. These and other definitions are equivalent when dealing with 
the Laplacian on the whole space 1^. The correct definitions for the operators defined on 
bounded domains admit several options that are being investigated at the moment, [19, 74]. 
The interest in these fractional operators has a long history in Probability and other applied 
sciences. The systematic study of the corresponding PDE models with fractional operators 
is relatively recent, and many of the results have been established in the last decade, see for 
instance the presentations [56, 66, 80, 82, 83], where further references can be found. 

A part of the current research concerns linear or quasilinear equations of elliptic type. 
This is a huge subject with well-known classical references. Even if it includes Free Boundary 
Problems, mainly of the Obstacle Type, see [27], it will not be discussed here in itself for lack 
of space, because we want to present free boundaries in evolution problems. 


Linear evolution processes of anomalous type 

We will study evolution models that arise as variants of the heat equation paradigm, includ¬ 
ing however long-range effects. Thus, the difference between the standard and the fractional 
Laplacian consists in taking into account long-range interactions (jumps) instead of the usual 
interaction driven by close neighbors. The change of model explains characteristic new fea¬ 
tures of great importance, like enhanced propagation with the appearance of fat tails at long 
distances: such tails great differ from the typical exponentially small tails of the standard 
diffusion, and even more from the compactly supported solutions of porous medium flows. 
Moreover, the space scale of the propagation of the distribution is not proportional to t 1 / 2 as 
in the Brownian motion, but to another power of time, that can be adjusted in the model; 
this is known as anomalous diffusion. Anomalous diffusion is nowadays intensively studied, 
both theoretically and experimentally since it explains a number of phenomena in several 
areas of physics, finance, biology, ecology, geophysics, and many others, which can be briefly 
summarized as having non-Brownian scaling, see e.g. [4, 51, 88]. The fractional Laplacian op¬ 
erators of the form (—A) CT / 2 , cr € (0, 2), are actually the infinitesimal generators of stable Levy 
processes [4, 10, 48]. The standard linear evolution equation involving fractional diffusion is 

£ + (-A»)=0, 

which is the main model for anomalous diffusion. The equation is solved with the aid of 
well-known tools, like Fourier transform. Posed in the whole space, it generates a semigroup 
of ordered contractions in L l (WL N ) and has the integral representation 


u 



K s (x 


z,t)f(z)dz, 
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where K s has Fourier transform K s {^,t) = e~^ 2st . This means that, for 0 < s < 1, the kernel 
K s has the self-similar form K s (x,t ) = t~ N ^ 2s F(\x\/t 1 ^ 2s ) for some profile function F = F s 
that is positive and decreasing, and it behaves at infinity like F(r) ~ r ~( N + 2s ) ; [18] . When 
s = 1/2, F is explicit: 

F 1/2 (r) = C (a 2 + r 2 )~^ N+l)/2 . 

If s = 1 the function K s =i is the Gaussian heat kernel, which has a negative square exponential 
tail, i.e., a completely different long-distance behavior. 

From the point of view of the present investigation, this fractional linear model and other 
variants, though useful and heavily studied, fail to meet the requirement that we are looking 
for, that is, generating free boundaries in the evolution. On the contrary, the solutions have 
large densities at long distances. 

3 Nonlinear and nonlocal diffusion models 

A main feature of current research in the area of PDEs is the interest in nonlinear equations 
and systems. There are a number of models of evolution equations that have been proposed as 
nonlinear counterparts of the linear fractional heat equation, and combine Laplace operators 
and nonlinearities in different ways. Let us mention the two popular in the recent PDE 
literature, and how they relate to our stated goal. 

• Model I. A quite natural option is to consider the equation 

dtu + (-Anu m )=0 (3.1) 

with 0 < s < 1 and m > 0. This is mathematically the simplest fractional version of the 
standard Porous Medium Equation (PME) dtu = A (u m ), that is recovered as the limit s —>• 1 
and has been extensively studied, cf. [5, 81]. We will call equation (3.1) the Fractional 
Porous Medium Equation, FPME, as proposed in our first works on the subject [53, 54], in 
collaboration with de A. de Pablo, F. Quiros and A. Rodriguez. From the point of view of our 
present study, this model is quite interesting because it offers a possible balance between a 
porous medium nonlinearity, which for m > 1 implies finite propagation and free boundaries, 
and a fractional-type diffusion where the speed on propagation is infinite and fat tails arise. 
The question is: which effect wins? In the papers [53, 54] we have answered the question 
in the sense of infinite propagation and no free boundaries (we are considering nonnegative 
solutions of the evolution problem). 

Nowadays, a reasonably complete theory of nonlinear diffusion with long-range effects has 
been developed for this model, see for instance [55, 84, 86, 20], or similar models like [50], but 
we have to look elsewhere for free boundaries. 

• Model II. A second model of nonlinear diffusion with fractional operators has been studied 
by the author in collaboration with Luis Caffarelli and does give rise to the occurrence of finite 
propagation and free boundaries. This alternative model is derived in a more classical way 
from the Porous Medium Equation since it is based on the usual Darcy law (i. e., the velocity 
of the particles is assumed to be the gradient of a pressure function, v = — Vp) with the 
novelty that the pressure is related to the density by an inverse fractional Laplacian operator, 
p = Ku. Putting these two facts into the continuity equation, dtu + V • (uv) = 0, the model 
takes the form 

ut = V(u V/Cu), (3-2) 
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a nonlinear fractional diffusion equation of porous medium type where /C is the Riesz operator 
that typically expresses the inverse to the fractional Laplacian, /C u = (—A )~ s u. This has been 
studied by Caffarelli and the author in [29, 30]. In dimension N = 1 this model was studied 
by Biler, Karch and Monneau [15, 14] as a model for the propagation of dislocations as 
proposed by Head [63]. Since it was extensively reported in the survey paper [82] we will 
only give point out to the established theory of weak solutions, and the lack of uniqueness 
results in several dimensions, and concentrate on the existence and properties of the free 
boundaries. Self-similar solutions exist and their existence and properties are very illustrative 
of the phenomenon of finite propagation. 

3.1 Finite propagation for Model II. Solutions with compact support 

One of the most important features of the porous medium equation and other related degen¬ 
erate parabolic equations is the property of finite propagation, whereby compactly supported 
initial data uq(x) gives rise to solutions u(x,t ) that have the same property for all positive 
times, i.e., the support of u(-, t) is contained in a ball B R(q(0) for all t > 0. One possible proof 
in the case of the PME is by constructing explicit weak solutions exhibiting that property 
(i.e., having a free boundary) and then using the comparison principle, that holds for that 
equation. Since we do not have such a general principle here, we have to devise a comparison 
method with a suitable family of “true supersolutions”, which are in fact some quite excessive 
supersolutions. The technique has to be adapted to the peculiar form of the integral kernels 
involved in operator K s . 

We begin with N = 1 for simplicity. We assume that our solution u(x, t) > 0 has bounded 
initial data uq(x) = u(x,to) < M with compact support and is such that uq is below the 
parabola a(x — b) 2 , a,b > 0, with graphs strictly separated. We may assume that uq is 
located under the left branch of the parabola. We take as comparison function 

U(x,t ) = a(Ct — (x — b)) 2 , 

which is a traveling wave moving to the right with speed C that will be taken big enough. 
Then we argue at the first point and time where u(x, t ) touches the left branch of the parabola 
U from below. The key point is that if C is large enough such contact cannot exist. The 
formal idea is to write the equation as 


u t = u x p x + up xx 


and observe that at the contact we have ut > Ut = 2 aC(Ct — x + b), while u x = U x = 
—2 a(Ct — x + 6), so the first can be made much bigger than the second by increasing C. The 
influence of p x and p xx as well as u is controlled, and then we conclude that the equation 
cannot hold if C is large enough. The argument can be translated for several dimensions. 
Here are the detailed results proved in [29]. 

Theorem 3.1. Let 0 < s < 1/2 and assume that u is a bounded solution of equation (3.2) 
with 0 < u(x, t) < L, and uq lies below a function of the form Uq{x) = Ae~ a \ x \ A, a > 0. If A 
is large then there is a constant C > 0 that depends only on ( N , s, a, L, A) such that for any 
T > 0 we will have the comparison 

u(x,t) < Ae ct ~ a ^ for all x € R N and all 0 < t < T. 
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If 1/2 < s < 1 a similar statement is true but C is not a constant but some increasing function 
of time. 

A similar finite propagation result is true in several space dimensions. The study of the 
corresponding free boundaries is an open topic. Let us mention in passing that much effort 
has been devoted in papers [29, 28, 31] to prove the basic sequence of regularity results, 
according to which initial data in L 1 give rise to solutions that belong to L p , 1 < p < oo for 
all positive time. Moreover, there exists a positive constant C such that for every t > 0 

sup \u(x,t)\ < C{N,s) t~ a ||itoll^i^ 

x£R n V 

with a = N/(N + 2 — 2s), 7 = (2 — 2 s)/((N + 2 — 2s). Finally, it is proved that bounded 
weak solutions u > 0 of the initial value problem are uniformly continuous on bounded sets 
of s < 1. Indeed, they are C a continuous with a uniform modulus. 

3.2 Self-similar solutions for Model II 

Next we want to study of the large time behaviour of this model following paper [30], since 
it shows some light on the way finite propagation works. The first step is constructing the 
self-similar solutions that will serve as attractors. Surprising result: their density will be 
compactly supported, while their pressure will be positive everywhere with the typical fat 
tails far away. 

Rescaling for the FPME. Inspired by the asymptotics of the standard porous medium 
equation, we define the rescaled (also called renormalized) flow through the transformation 

u(x, t) = (t + l)~ a v(x/(t + 1)^, r) (3.3) 

with new time r = log(l+t). We also put y = x/(t + 1)P as rescaled space variable. In order to 
cancel the factors including t explicitly, we get the condition on the exponents a + (2 — 2 s)/3 = 
1. Here we use the homogeneity of K, in the form (Ku)(x,t) = t~ a+2s ^ {Kv){y, r). From 
physical considerations we also impose the law that states conservation of (finite) mass, which 
amounts to the condition a = N/3, and In this way we arrive at the precise value for the 
exponents: 

P = l/(JV + 2 —2s), a = N/{N + 2-2s). 

We also arrive at the nonlinear, nonlocal Fokker-Planck equation 

v T = \7y(v(\7 y lC(v) + py)) (3.4) 

which is the equation for the renormalized Bow. In all the above calculations the factor (t +1) 
can be replaced by t + to for any > 0 , or even by plain t. 

Stationary renormalized solutions. It is important to concentrate on the stationary states 
of the new equation, i. e., on the solutions V(y) of 

V y ■ (W y (P + a\y\ 2 )) = 0, with P = K{V). 

where a = (3/ 2 , and /3 is defined just above. Since we are looking for asymptotic profiles of the 
standard solutions of the FPME we also want V > 0 and integrable. The simplest possibility 
is integrating once to get 

V S7 y (P + a\y\ 2 )) = 0 , P = K(y), V > 0 . ( 3 . 5 ) 
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The first equation gives an alternative choice that reminds us of the complementary formula¬ 
tion of the obstacle problems. Indeed, if we solve the Obstacle problem with fractional 
Laplacian we will obtain a unique solution P(y) of the problem: 

P>f, V = (-AyP>0 ; 
either P = f or V = 0. 

with 0 < s < 1. In order for solutions of (3.6) to be also solutions of (3.5) we have to choose 
as obstacle f(y) = C — a\y\ 2 , where C is any positive constant and a = (3/2. Note that 
—A/ = 2 Na = a. For uniqueness we also need the condition P —> 0 as \y\ — > oo. Fortunately, 
the corresponding theory had been developed by Caffarelli and collaborators, cf. [26, 75]. The 
solution is unique and belongs to the space H~ s with pressure in H s . Moreover, it is shown 
that the solutions have P £ C 1,s and V € C 1-5 . 

Note that for C < 0 the solution is trivial, P = 0, V = 0, hence we choose C > 0. We 
also note the pressure is defined but for a constant, so that we could maybe take as pressure 
P = P — C instead of P so that P = 0; but this does not simplify things since P — > 0 implies 
that P —» — C as |y| — > oo. Keeping thus the original proposal, we get a one parameter family 
of stationary profiles that we denote Vc(y)■ These solutions of the obstacle problem produce 
correct weak solutions of the fractional PME equation with initial data a multiple of the Dirac 
delta for the density, in the form 

U c (x,t) =t- a V c (\x\t-P). 

It is what we can call the source-type or Barenblatt solution for this problem, which is a 
profile V > 0. It is positive in the contact set of the obstacle problem, which has the form 
C = {|y| < R(C)}, and is zero outside, hence it has compact support. 

On the other hand, the rescaled pressure P(\y\) is always positive and decays to zero as 
|y| —> oo according to fractional potential theory, cf. Stein [79]. The rate of decay of P as 
|y| —> oo turns out to be P = 0(\y\ 2s ~ N ), a fat tail. 

Exact calculation of density profiles. Biler, Karch and Monneau [15] studied the existence 
and stability of self-similar solutions in one space dimension. Recently, Biler, Imbert and 
Karch [14] obtain the explicit formula for a multi-dimensional self-similar solution in the form 

U(x,t) = cR- Q (l - xH~^l N f~ s 

with a = N / (N + 2 — 2s) as before. The derivation uses an important identity for fractional 
Laplacians which is found in Getoor [62]: (—A)°'/ 2 (l — y 2 ) a + 2 = if a £ (0,2], Here we 

must take a = 2(1 — s). According to our previous calculations A P = —a on the coincidence 
set, hence c\ = a/K a ^. Let us work a bit more: using a similar scaling to (3.3), we arrive at 
the following one-parameter family of self-similar solutions 

U(x,t] Ci) = t“"(Ci - h x 2 t~ 2a i N )\- s 

where k\ = ^ and Ci > 0 is a free parameter that can be fixed in terms of the mass 

of the solution M = f U(x, t; C\) dx. This is the family of densities that corresponds to the 
pressures obtained above as solutions of the obstacle problem. Let us finally mention that 
explicit formulas for self-similar solutions to related fractional equations have also recently 
been obtained in [64], 
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3.3 Asymptotic behaviour for Model II 

The next step is to prove that these profiles are attractors for the rescaled flow. The studay 
is done in terms of the rescaled flow that is much more workable than the original evolution. 
The following result is proved in [30]. 

Theorem 3.2. Let u(x,t ) > 0 be a weak solution of the Cauchy Problem for equation (3.2) 
with bounded and integrable initial data such that uq > 0 has finite entropy, i.e., 

I (v K.(v) + (3\y\ 2 v) dy < oo . 

J R N 

Let v(y,r) be the corresponding rescaled solution to (3.4). As r —>• oo we have 

v(-,t) Vc(y) in L 1 (M. N ) and also in L°°(R N ). 

The constant C is determined by the rule of mass equality: f Rjv v(y , r) dy = J Rjv Vc(y) dy. In 
terms of function u, this translates into the convergence as t —>• oo 

u(x,t) — Uc(x,t) 0 in L 1 (R N ), t a \u(x, t) — Uc(x, t)| —S > 0 uniformly in x. 


The constants a and /3 are the self-similar exponents defined before. The proof uses entropy 
dissipation methods that adapt perfectly to the problem. Sharper results on this convergence 
are given in [43] where more references can be found. See also below in this respect. 
Extension of the nonlinear fractional models. 1) There is way of generalizing the two 
previous Models I and II, to accept two exponents m and p: 

d t u + V(u"* _ 1 V(- A)~ s u p ) = 0 , 

so that the comparison of both models happens on symmetric terms. The question we want 
to solve is deciding between finite and infinite propagation in terms of the exponents m and 
p. Progress on this issue is reported in [77, 78]. 

2) Limits. The limit of Model II when s —> 1 is quite interesting since one obtains a variant 
of the equation for the evolution of vortices in superconductivity studied in [47], E [58, 69, 2]. 
The understanding of this limit has been done in collaboration with Sylvia Serfaty [73], and 
is related to work by Bertozzi et al. on aggregation models [11, 12]. This is a good point to 
connect to the next issue. On the other hand, the limit m —> oo of Model 1 leads to a free 
boundary problem of the Mesa Type reported in [84, 85]. 


4 Aggregation and swarming 

The collective behavior of individuals in animal grouping, also called swarming, has recently 
been modelled by mean field PDEs, see [67] and the references therein. They consist in 
continuum descriptions of Individual Based Models (IBMs), see [3, 72, 65] for applications 
to fish schools and software simulation of swarms. The motion of individuals in IBMS is 
governed by systems of ODEs imposing an asymptotic speed, cruise speed, for particles. One 
of the most celebrated models of this kind was introduced in [57]. The authors propose to 
introduce a social force between individuals based on the empirical observation that there are 


two basic mechanisms of interaction: a inner repulsion zone to keep a comfort area or to avoid 
collisions and an outer attraction area since individuals do want to socialize and be close to 
each other. These two tendencies are modeled via a pairwise potential directed along the 
position between the particles and depending on the interparticle distance. Therefore, this 
model consists in 


dxi 

dt 


= Vi 


dvj 

dt 


= avi - (5vi\vi 


^2 VlT(xi - Xj) , 
j¥=i 


(4.1) 


where Xi,Vi € M. N ,i = 1 ,M are the positions and velocities of the individual particles and 
a, (3 are effective values for self-propulsion and friction forces, see [57] for more discussions. 
These models lead to kinetic equations in the limit of large number of particles as derived in 
[40]. Some of the patterns that the simulation of these IBMs produce are the so-called flock 
patterns. The flock solution corresponds to the case in which all particles move with the same 
velocity vector Vi = vo with |uo| 2 = a//3 while the positions x* € M. N , i = 1 ,,M satisfy that 
the forces are balanced, that is, 


T: VIE (Xi — Xj) = 0 . 


The flocking solution corresponds then to the translational movement of the shape formed by 
the positions Xi £ M. N ,i = 1 with velocity vector vq. When the number of particles 

tends to infinity, these discrete flocking solutions are approximating continuum density profiles 
p satisfying VIE * p = 0 at least over the support of p. In fact, the connection between the 
second-order microscopic model (4.1) and the first-order model given by 


dxi 

dt 


Ev^.-^), 


(4.2) 


has been clarified in [1, 42]. There, it is shown that the stability of flock solutions for (4.1) 
is equivalent to the stability of the spatial shape of the flock as steady solution to (4.2). A 
similar result at the level of the continuum equations is lacking. However, it makes sense 
to find all the spatial configurations leading to flocks or equivalently all stationary solutions 
to (4.2). Among those, the most stable configurations are of special interest, in view of the 
connection through the stability of the solutions at least at the discrete level. The system 
(4.2) has a natural energy that is dissipated along its evolution given by 

E[X!, ... ,x M \ = 2 ^2 W ( x i~ x j)- 

The continuum version of this energy, suitably scaled, leads to 

E\fA = \ f W(x — y)dp(x)dp(y) (4.3) 

z Jr n xM. n 

where p is a probability measure. This energy is dissipated along the evolution of the contin¬ 
uum mean-field limit of the system of ODEs in (4.2), called the aggregation equation 


d t p = div(p(VlT * p )), 


(4.4) 
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where p is the density function of the measure p. Therefore, minimizers of the energy (4.3) 
on the set of probability measures should be among the most stable configurations. Finally, 
let us comment that both the nonlocal diffusion equation (3.2) and its rescaled version (3.4) 
are particular instances of the aggregation equation (4.4) for suitably chosen potentials W. 
We will see in the rest of this section a summary of the results concerning the minimizers of 
the interaction energy (4.3), specially in its relation to free boundary problems, particularly 
to the obstacle problem. 

Connection to optimal transport. We start by noting that the discrete system (4.2) has 
the structure of a gradient flow with respect to the Euclidean distance. In fact, this structure 
is shared by the continuum equation (4.4) in a suitable sense. We need to be precise about 
the definition of distance between probability measures. We recall that P(R N ) is the set of 
Borel probability measures on R^ and we denote by B(R N ) the family of Borel subsets of 
R N . The support of a measure p £ P(R N ) is the closed set defined by 

supp (p) := {x £ R N : p(B(x , e)) > 0 for all e > 0} . 

A family of distances between probability measures has been classically introduced by means 
of optimal transport theory, we will review briefly some of these concepts, we refer to [87] 
for further details. A probability measure it on the product space R^ x R^ is said to be a 
transference plan between p £ P(R N ) and v £ P(R N ) if 

7 t(A x M n ) = p(A) and ^(R^ x A) = is (A) (4.5) 

for all A £ B(R n ). If p, is £ P(U N ), then 

I l(p, u) := {n £ P(R n x R w ) : (4.5) holds for all A £ ^(R^)} 

denotes the set of admissible transference plans between p and is. Informally, if 7 r £ II(p, is) 
then dn(x, y) measures the amount of mass transferred from location x to location y. We 
recall that the Euclidean Wasserstein distance d ,2 between two measures p and is is defined by 


dj(p,is)= inf 

7ren(/i,v) 


\x - y\ 2 dir(x,y) 


Note that d 2 {p, is) < oo for p,is £ P 2 (R N ) the set of probability measures with finite moments 
of order 2. It is classical that the equation (4.4) is a gradient flow of the energy (4.3) in 
the space of measures 7^2 (R' V ) endowed with the distance d 2 , see [44, 45]. This property 
is essential to obtain many qualitative properties of the solutions and to deal with possible 
singular interaction potentials, see [11, 38, 39, 9, 12, 6, 7, 8, 34] and the references therein. 

In our context of minimizing the energy functional (4.3), this means that we look for local 
minimizers in cfo- (Local) Minimizers of (4.3) should correspond to equilibrium configurations 
for the evolution equation obtained by steepest descent of the energy. However, being a func¬ 
tional on probability measures, the steepest descent has to be understood in the Wasserstein 
sense as in [70] by writing (4.4) as 


dpt 

dt 


div 



div (p t Vipt) 


x £ R^ , t > 0. 


(4.6) 
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This evolution equation can make sense in the set of probability measures depending on 
the regularity of the potential W. However, if the potential is singular at the origin the 
well-posedness only happens in L p spaces, see [13, 38, 36] for more details. For steady config¬ 
urations, we expect S7ip = 0 on the support of p, due to the formal energy dissipation identity 
for solutions, i.e., 

4^N = - [ IV-0*1 2 dpt 

dt J r n 

where pt is any solution at time t of (4.6) and ipt = W * pt its associated potential. Therefore, 
the points on the support of a local minimizer p of the energy E should correspond to critical 
points of its associated potential ip. This fact is made rigorous in: 

Proposition 4.1 ([7, Theorem 4]). Assume that W is a non-negative lower semi-continuous 
function in Lj oc ( l w ). If p is a d 2 -local minimizer of the energy, then the potential ip satisfies 
the Euler-Lagrange conditions given by 

(i) ip(x) = (W * p){x) = 2E\n] p-a.e. 

(ii) ip(x) = {W * p){x) < 2 E\p\ for all x £ supp(p). 

(in) ip(x) = (W * p)(x) > 2 E\p] for a.e. x £ M. N . 

These conditions simplify to 

ip{x) = (W * p){x) = 2 E[p] for a.e. x £ supp(^), 

ip{x) = (W * p)(x) > 2 E[p] for a.e. x £ \ supp(/i), (4.7) 

if /i is absolutely continuous with respect to the Lebesgue measure. These Euler-Lagrange 
conditions can also be interpreted in terms of game theory, see [33, 16, 52], In fact, they imply 
that the potential ip achieves its global minima on the support of p,, and therefore p is also 
characterized as 

/ ipix)dp(x) = min / ip(x)di / (x). 

Jrn vGV2(R n )Jrn 

In game theory the probability measure p is interpreted as the strategy chosen by the agent 
x , the potential is interpreted as encoding the interactions between the agents, and the mini¬ 
mization is thought about the optimal condition (the game) that the agents are looking for. 
Therefore, the Euler-Lagrange conditions are interpreted in game theory as the lack of interest 
of any player in the game to change of strategy leading to the concept of Nash equilibria. As 
a summary, local minimizers of the potential energy are among Nash equilibria within the 
game theoretical viewpoint. 

Let us concentrate on the quest of local minimizers for repulsive/attractive potentials of 
the form 

I rp I (L I 'T* I ^ 

W{x) = u. - LL , (4.8) 

a b 

with a > b > — N to have a repulsive at the origin and attractive at infinity potential which 
is locally integrable. We allow a or b to be zero with the understanding that |x|°/0 ~ log |x|. 

As mentioned before, the case a = 2 and b = 2 — N, black point in Figure 1, is of particular 
interest. It corresponds to Newtonian repulsion and it has been repetitively rediscovered that 
the unique, up to translations, global minimizer of the interaction energy is the characteristic 
function of a Euclidean ball with suitable radius [61, 71, 46]. The uniqueness, up to transla¬ 
tions, of the global minimizer for more singular than Newtonian repulsion, a = 2, b = 2s — N 
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Figure 1: Parameters for minimizers of the interaction energy with potential of the form (4.8) 
given by densities, see the text for a full explanation. 

with 0 < s < 1 (red dashed line in Figure 1), was obtained in [30] via the connection to a 
classical obstacle problem, see also [46]. This strategy was also used in [73] to treat again the 
case s = 1 for the evolution problem as in [12]. If the parameters are in the range of the green 
line in Figure 1, i.e. a > 2 or 2 — N < a < 2 with b = 2 — N, the existence and uniqueness 
of compactly supported radial minimizers of the interaction energy was obtained in [59, 60]. 
Summing up, the solution support is limited by a free boundary, whose location (a sphere) is 
to be determined as part of the problem. 

Moreover, the main result in [37] shows that in the range a > 0 with b = 2s — N with 
0 < s < 1, grey area in Figure 1, local minimizers in c ?2 of the interaction energy (4.3) are 
absolutely continuous with respect to the Lebesgue measure, and their density function lies 
in (M. n ) when s = 1 and in (M ;¥ ) when s € (0,1). These results are obtained again 
by exploiting the connection between the Euler-Lagrange conditions for local minimizers and 
classical obstacle problems. They can be generalized for repulsive at the origin potentials 
behaving like \x\ b with b = 2s — N, 0 < s < 1, and under quite general conditions of the 
behavior of the potential W outside the origin. In all these cases, we know that the local 
minimizers are compactly supported and determined by a free boundary whose regularity is 
still unknown in full generality. 

Let us conclude this part by mentioning that whenever we can fully characterize the 
unique, up to translations, global minimizers of (4.3), they are the long-time asymptotics for 
the aggregation equation (4.4), see [12, 29, 30, 73] for the known different cases (and Section 3 
above). These cases are essentially reduced to a = 2 and b = 2s — N with 0 < s < 1, and even 
decay rates are known in one dimension [41, 43]. Finally, let us remark that finding sufficient 
conditions on the interaction potential W to ensure the existence of compactly supported 
global minimizers is a very interesting question that has been recently solved in [32], see also 
[76, 49, 35]. However, we do not know yet about their properties in terms of asymptotic 
stability for the aggregation equation (4.4). 

More on the connection to free boundary problems and the obstacle problem. 

We will now focus on a topic already introduced in Section 3 (b) above. A classical obstacle 
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problem [26] for a differential operator C consists in finding a function : M. N —»• R such that 


($-/)(£$-<?) = 0 

4> > / on R^, (4.9) 


where /, g : R w —>■ R are just given functions. Here, / is the obstacle to <h, and thus whenever 
/ does not touch the obstacle, then the first condition implies that <f> should be a solution to 
the PDE: = g. To make the connection to the Euler-Lagrange condition clearer, let us 

take the particular case of the classical potential encountered in semiconductors [25, 24] given 
by 


W (x) = cn 


x 


2—TV 


2 - N 



with N > 1 . Here, cjsr is the constant such that the first part of the potential is the fundamental 
solution of the Laplacian operator A in IV > 1. Assume that p £ L 1 n L°°(R N ) is a cfo-local 
minimizer of the energy associated to this potential, then the potential associated to p, given 
by '0 = W * p is well defined as a function that satisfies in the sense of distributions 


Aip = AW * p = p + N, 

and thus, Aifr > N in R^ and Aifi = N outside the support of p. In conclusion, the potential 
ip associated to the local minimizer satisfies a classical obstacle problem in R^ with operator 
C = A, g = IV, and / = 2 E[p\ due to (4.7). Let us also remark that the obstacle problem 
(3.6) associated to the scaled version of the fractional diffusion equation (3.4) can also be 
recast in the form (4.9). We will do this next but in more generality. 

The basic idea for the Laplacian can be generalized and used to get properties on the 
potentials associated to local minimizers of a larger family of interaction potentials W without 
resorting to the local minimizers themselves. Let us review quickly the main strategy in [37]. 
The key fact is that if the repulsive strength at zero of the potential is given by \x\ b /b, 
b = 2s — N with 0 < s < 1, then one can show that the potential if; = W * p associated to 
any c^-local minimizer p is a continuous function, see [37]. This allows us to deduce that for 
any point xq € supp (p) 

ij}{x) > V’(tco) for all x £ B £ (xq) (4-10) 

holds. Here, e is the size of the ball in cfo distance where p is a local minimizer. Moreover, 
for the points in the support, we also get 


ip(x) = ipixo) in B £ (xo) n supp(^). (4.11) 

Since ^gj- in this range is the fundamental solution, up to constants, of the fractional diffusion 
operator C = — (—A) s , we can write in the distributional sense that 

{-A) s ^ = p + (-A) s ^ *p inV'(R N ). 

a 

In particular, since p is a non-negative measure, we deduce 

(-A ) s ip > (-A) s — * p in B e (x 0 ). 
a 
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Furthermore, if x £ B e (x o) is such that ip(x) > if{xo), (4.11) implies that x ^ supp(p). We 
deduce 

(-A) s -i l> = (-A) s -^-L * p in T>'(R n ) in B e (x 0 ) n {if > if{x 0 )} . 

Collecting (4.10)-(4.11), we have, at least formally, the following result: 

Theorem 4.2 ([37]). For all xq £ supp(p), the potential function if associated to a d 2 -local 
minimizer of the energy (4.3) is equal, in B e (x o), to the unique solution of the obstacle problem 

<P > Co, 

(-A)V > ~F{x), 

(-A)V =~F{x), 

V = Vh 

where Co = if(xo) and F(x) = —(—A) s -^|- * p. 
mizer of the energy p as p = —A if + F. 

Most of the results mentioned above follow this strategy to make use of the regularity and 
existence theorems for obstacle problems available in the literature, see [22, 23, 17, 75, 26]. 
Let us remark that in the case of a power law potential (4.8) with a = 2 and 0 < s < 1, F[x) 
is constant. The range corresponding to 0 < s < 1 and a = 2 corresponds to the obstacle 
problem treated by Caffarelli and Vazquez in [30]. As before, Theorem 4.2 is true for more 
general potentials behaving as the repulsive power |x| fc at the origin with suitable growth 
conditions at infinity, see [37]. 

It is an open question how to use obstacle problem techniques in cases in which the 
repulsion is not too singular, for instance for power-law potentials when (4.8) with a = 2 and 
b = 2s — N with now s > 1. The main problem we face is that the inverse operator of the 
convolution is no longer as nice as it was for fractional diffusion. However, we are currently 
investigating the case when s = 1 + e with 0 < e < 1, i.e., in the range corresponding to the 
blue line in Figure 1. In this particular case, we consider an intermediate obstacle problem 
by defining fk to be the solution to 

*>ci, V = (—A) e \k > 0; 

either 'F = ci or V = 0 , 

which again has nice regularity properties since 0 < e < 1. Now, we define our potential if 
as the Newtonian potential associated to 'F, that is, we solve \F = —Aip. The idea is that 
by integrating twice the solution of the obstacle problem (4.13), we get a candidate to be a 
solution to (4.12). 


in B e (x 0 ) 
in B e (x 0 ) 

in B e (x o) n {ip > Co} 
on dB £ (x 0 ), 


(4.12) 


Furthermore, we can recover the local mini- 
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